Unexpected Density Fluctuations in Jammed Disordered Sphere 

Packings 

Aleksandar Donev/'^ Frank H. Stillinger,^ and Salvatore Torquato^'^'^'Q 

^Program in Applied and Computational Mathematics, 
Princeton University, Princeton NJ 08544 
^PRISM, Princeton University, Princeton NJ 08544 
■^Department of Chemistry, Princeton University, Princeton NJ 08544 

Abstract 

We computationally study jammed disordered hard-sphere packings as large as a million particles. 
We show that the packings are saturated and hyperuniform, i.e., that local density fluctuations 
grow only as a logarithmically-augmented surface area rather than the volume of the window. The 
structure factor shows an unusual non-analytic linear dependence near the origin, S{k) ~ \k\. In 
addition to exponentially damped oscillations seen in liquids, this implies a weak power-law tail in 
the total correlation function, h{r) ~ —r~'^, and a long-ranged direct correlation function. 
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The characterization of local density fluctuations in many-particle systems is a problem 
of great fundamental interest in the study of condensed matter, including atomic, molecular 
and granular materials. In particular, long-wavelength density fluctuations are important to 
such diverse fields as thermodynamics, granular flow, and even cosmology (see Ref. 3| and 
references therein). Previous work by some of us jl'j was concerned with the quantitative 
characterization of density fluctuations in point patterns, and in particular, those in which 
infinite wavelength fluctuations are completely suppressed, i.e., the structure factor S(k) 
vanishes at the origin. In these so-called hyperuniform (or superhomogeneous P]) systems, 
the variance in the number of points inside a large window grows slower than the volume of 
the window, typically like the window surface area. Most known examples of hyperuniform 
systems are either ordered lattices or quasi- crystals . An important open problem is 
the identification of statistically homogeneous and isotropic atomic systems (e.g., glasses) 
that are hyperuniform. 

For equilibrium liquids and crystals, S{k = 0) is proportional to the isothermal compress- 
.biUty and i= thu= po.jUve. SMctiy ,a^^e, sphere packings Q a.e rigo.ou^ly _s.ible 
(and non-shearable) j^, but they are also nonequilibrium systems. In Ref. Pl it was con- 
jectured that all saturated [5] strictly jammed packings are hyperuniform. Of particular 
importance to understanding both glasses and granular materials are disordered jammed 
packings, and in particular the maximally random jammed (MRJ) state j^, which is the 
most disordered among all strictly jammed packings [3]. The MRJ state for hard-particle 
packings is related to the view of jamming as a rigidity transition and/or dynamic arrest in 
both granular 0] and glassy materials 0]. Previous studies have identified several different 
diverging length scales at the rigidity jamming transition for systems of soft spheres (see Ref. 
0] and references therein), indicating a kind of second-order phase transition at the jamming 
point. Hyperuniformity involves an "inverted critical phenomenon" in which the range of the 
direct correlation function c(r) diverges It is therefore of great interest to test whether 
disordered jammed sphere packings are hyperuniform. In this Letter, we demonstrate that 
MRJ packings are indeed hyperuniform and saturated. Moreover, we observe an unusual 
non-analytic structure factor S{k) ~ \k\, or equivalently, a quasi-long ranged negative tail 
of the total pair correlation function h{r) ~ — r~^, just as found in the early Universe P|. 

We prepare jammed packings of hard spheres under periodic boundary conditions using 
a modified Lubachevsky-Stillinger (LS) packing algorithm [l^, as detailed in Ref. jll| . 
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The generated disordered sphere packings typically have volume fractions in the range = 
0.64 — 0.65, and to a good approximation the packings should be representative of the MRJ 
state. For this study, we have generated a dozen packings of = 10^ and = 10^ particles 
jammed up to a reduced pressure of 10^^ using an expansion rate of 10~^ with cf) ^ 0.644. 
Generating such unprecedented one-million-particle packings was neccesary in order to study 
large-scale density fluctuations without relying on dubious extrapolations. 

The packings generated by the LS and other algorithms have a significant fraction 
(~ 2.5%) of rattling particles which are not truly jammed but can rattle inside a small 
cage formed by their jammed neighbors These rattlers make a negligible contribution 
to the mechanical properties of the system, including the pressure, and can be removed suf- 
ficiently close to the jamming point. However, they are important when considering density 
fluctuations. Removing the rattlers will produce small but observable long-wavelength den- 
sity fluctuations. Assuming that the rattlers are more or less randomly distributed among 
all particles, a hyperuniform packing from which the rattlers are then removed would have 
S'(O) ~ 0.025 > 0. Similarly, the hyperuniformity could be destroyed by randomly filling 
large-enough voids with additional rattlers. It is therefore important to verify that the 
jammed packings are saturated, i.e., that there are no voids large enough to insert addi- 
tional rattlers. Figure ^ shows the complementary cumulative pore-size distribution 12] 
which gives the probability that a sphere of diameter 5 could be inserted into the 
void space, with and without the rattlers. Clearly there is no room to insert any additional 
rattlers; the largest observed voids are around 5max ~ 0.8/^. The algorithm used to produce 
the packings appears to fill all void cages with particles. 

It is very difficult to study long-wavelength density fluctuations accurately in 3D com- 
puter simulations. When periodic boundary conditions apply with a periodic box of length L, 
particle correlations can only be studied up to a distance L/2, and there are large finite-size 
corrections for distances comparable to L. Additionally, as we show later, strong statis- 
tical fluctuations appear due to finite system size, making it necessary to use even larger 
systems to measure pair correlations at large distances. In reciprocal space, S{k) can only 
be measured for k > 27r/L, with large discretization errors for the smallest wavevectors. 
To overcome these finite-size effects, it was necessary to generate a packing of one million 
particles. 

Consider a large isotropic three-dimensional packing of A^ hard spheres of diameter D, 
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Figure 1: (Color online) The cumulative pore-size distribution F{5) for a (single) packing with 
N = 10^ par ticles, with and without the rattlers. The method of trial spheres with 2 • 10^ trials 
was used 12]. A very similar F{6) with cutoff around 5max ^ 0.8D is observed when N = 10^, 
when rattlers are present. The cutoff is however not as sharply defined as it is for the FCC crystal, 
shown for comparison, since the exact size of the largest available void (cavity) fluctuates between 
realizations. 

with average number density p = N/V and average volume fraction (j) = TrpD^/Q. We 
employ the usual pair correlation function g2{x = r/D) or the total correlation function 
h{x) = g2{x) — l in real space, or the equivalent Fourier representation given by the structure 
factor 

S{K = kD) = l + 240 r ^-^^^^x^h{x)dx. 

Jo Kx 

Of particular interest are the moments of /i(x), (x") = x'^h{x)dx. Computer-generated 
packings are always finite, and thus binning techniques must be used to obtain probability 
densities like h. Accordingly, we prefer to use the more readily measurable excess coordina- 
tion 

AZ (a;) = 1 + 240 f w^h{w)dw. 
Jo 

This is the average excess number of points inside a spherical window of radius xD centered 

4 



at a particle, compared to the ideal-gas expectation 80x^. Any integral containing h{x) can 
easily be represented in terms of AZ {x) using integration by parts. For the structure factor 
we get S (K) = lim^j^oo S{K, R), where 

This has quadratic behavior near k = when expanded in a Taylor series, 

poo 

S (K) ^ 8(0) + ^ / X [AZ{x) - 5(0)] dx, (2) 
-J Jo 

where S{0) = AZ (x — oo) vanishes for a hyperuniform system. For large x, an explicit 
finite-size correction of order needs to be applied to the infinite-system excess coor- 
dination, AZ{x) ~ S{0) [1 — 8(f)x^/N] |l3^, as it is clear that the excess coordination must 
vanish for windows as large as the whole system. 

Figure 121 shows S{k) for the simulated packings as obtained via a direct Fourier transform 



(DFT) of the particle positions, 5'(k) = N ^ J2iLi ^^P (^^ ■ r 



where k is a reciprocal 



lattice vector for the periodic unit cell . To obtain an approximation to the radially sym- 
metric infinite-system S{k), we average over the reciprocal lattice vectors inside a spherical 
shell of thickness 2tt/L. Using Eq. together with a numerical (truncated) AZ{x) quickly 
gives S{k) over a wide range of wavelengths. However, the behavior near the origin is not 
reliable since it depends on the analytic extension for the tail of AZ{x). The results of the 
DFT calculations are shown in Fig. |21 and they closely match the one obtained from AZ{x) 
for wavelengths smaller than about 20 diameters. 

F„.e0de.on...« .Ka. .Ke ..u.a« pacM,, . .deed ,„^o^ R .o w.Mn 
S{0) < 10~^, as conjectured in Ref. The behavior of S{k) near the origin is very 

surprising. The observed S{k) follows closely a non- analytic linear relationship well- 
fitted by S{K) ^ 6.1 • 10~^ + 3.4 ■ IQ-'^K over the whole range K/27r < 0.4. By contrast, 
analytic quadratic behavior is observed for a liquid sample at = 0.49, as shown in the 
figure. Theoretical finite-size corrections to the small-fc behavior of S(k) have only been 
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considered for relatively low-density liquid systems with relatively small N and do 
not appear useful for our purposes. Although estimating the corrections to the DFT data 
analytically is certainly desirable, such corrections appear to be rather small at least for the 
well- under stood liquid at = 0.49. Comparison among the different = 10^ samples shows 
that statistical fluctuations in S{k) near the origin are very small. 



Equation Q shows that if h is truly short-ranged, the structure factor must be analytic 
(i.e., an even power of k, usually quadratic), near the origin. Our numerical observations 
point strongly to a linear S{K) for small K. It is interesting to note that such non-analytic 
behavior is assumed in the so-called Harrison- Zeldovich power spectrum of the density fluc- 
tuations in the early Universe and has been verified experimentally with high accuracy. 
If this observation S{K) ~ \K\ survives simulations of even larger jammed hard sphere 
systems, using a variety of packing algorithms, it would imply a negative algebraic power- 
law tail h{x) ~ —x~^ uncharacteristic of liquid states. Such a quasi-long range correlated 
decrease in the density around a particle is typically only seen in systems with long-range 
interactions. A long-ranged tail must appear in the direct correlation function c(r) for a 
strictly hyperuniform system due to the divergence of c(0), in a kind of "inverted critical 
phenomenon" Q]. Such a tail is uncharacteristic of liquids where the range of c(r) is sub- 
stantially limited to the range of the interaction potential. The direct correlation function 
can numerically be obtained from its Fourier transform via the Ornstein-Zernike (OZ) equa- 
tion, c{k) = (n/Gcf)) [S{k) — 1]/S{k), and we have shown it in the inset in Fig. 121 along 
with the corresponding Percus-Yevick (PY) anzatz 0| for c(r) at = 0.49 which makes 
the approximation that c(r) vanishes outside the core. Two unusual features relative to the 
liquid are observed for our jammed packing. First, there is a positive (5-function at contact 
corresponding to the Z = 6 average touching neighbours around each jammed particle 
Second, there is a relatively long tail outside the core, the exact form of which depends on 
the behavior of S{k) around the origin [iflj ]. 

The numerical coefficient in the power-law tail in h{x) is very small, AZ{x) ~ 4.4-10~^x~^, 
and cannot be directly observed, as we will show shortly. It is however possible to observe its 
effect on large-scale density fluctuations. For monodisperse hard sphere systems it suffices 
to focus only on the positions of the sphere centers and consider density fluctuations in point 
patterns. Following Ref. consider moving a spherical window of radius R = XD through 
a point pattern and recording the number of points inside the window N{X). The number 
variance is exactly ll, 

a'{X) = (N'iX)) - {N{X)f 

= ^[{2XfAZo{2X)-AZ2{2X)] 

where AZn{x) = w"'AZ{w)dw denotes a running moment of AZ. Asymptotically, for 
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large windows, in an infinite system with analytic S{k), cr^(X) ^ AX^ + BX^, where A = 
8(f) {1 + 2A(j){x'^)) = 8(j)S{0) is the volume fluctuation coefficient, and B = — 1440^ (a;^) = 
60AZo (x — > oo) is the surface fluctuation coefficient. When a non-integrable power-law tail 
exists in AZ{x), asymptotically the "surface" fluctuation coefficient contains an additional 
logarithmic term, B{X) = Bq + C \n X . Such a logarithmic correction does not appear for 
any of the examples studied in Ref. Explicit finite-size effects for non-hyperuniform 
systems yield a correction A{X) = 8(j)S{0) [l~8(j)X^/N] [2^. Figure El shows numerical 
results for the number variance as a function of window size, along with the predicted 
asymptotic dependence, including both the logarithmic and A^~^ corrections 21j. Both 
corrections need to be included in order to observe this close a match between the data and 
theory. The constants >S'(0) and C were obtained from the linear fit to Sjk), while Bq ^ 1.02 
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was obtained by numerically integrating AZ( , clS explained shortly 

We now turn our attention to real space to observe directly the large-distance behavior 
of h or equivalently AZ. For equilibrium liquids with short-ranged potentials it is expected 
that the asymptotic behavior of h{x) is exponentially damped oscillatory j^, 2d\, of the 
form 

C 

h{x) ~ — exp (— x/^) cos [Ko{x — Xq)] . (3) 

However, it is not clear whether the decay is still exponential for glass-like nonequilibrum 
jammed systems. Previous studies have looked at much smaller systems, where explicit 
finite-size effects dominate, and also focused on the liquid phase Il3| . Figure lU shows the 
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numerical AZ{x) along with a relatively good exponentially damped oscillatory fit |2^ 
AZix) ^ 5.47x exp(— a;/1.83) cos(7.56a; — 2.86) over the range 5 < x < 15. It would be 
desirable to look at larger x and, in particular, directly observe the long-range inverse power 
tail predicted from the linear behavior of Siji). The use of cubic periodic boundary conditions 



implies that pair distances up to Xmax = \/T^N/2A(j) 50 can be studied. However, it is 
important to point out that it is not possible to measure the pair correlations for x > 15 
due to statistical variations among finite systems, estimated to lead to an uncertainty of the 
order 5Z(x) ~ a{x)/^/N. In fact, within the range of validity of the observed AZ(x) the 
damped oscillatory fit is perfectly appropriate. We smoothly combined the actual numerical 
data for x < 10 with the fitted decaying tail for x > 10, and numerical integration of this 
smoothed AZ{x) gives Bq ^ 1.02 ±0.02, as used in producing Fig. El This smoothed AZ{x) 
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was used to obtain S{k) via Eq. (jT)) when producing Fig. |21 

We have given computational results for a million-particle jammed disordered hard sphere 
packing demonstrating that it is nearly hyperuniform and saturated. However, there are 
many open fascinating questions. Can a geometrical significance be attached to the period 
of oscillations Kq in the jamming limit 26|, or to the cutoff of F{6)7 We believe that 
the strict jamming and saturation conditions demand hyperuniformity of our packings. We 
conjecture that the observed non-analytic behavior of S{k) ~ is a direct consequence 
of the condition of maximal disorder on the jammed packing. The exponent p appears to 
increase with increasing order: It approaches infinity for ordered lattices, is two for perturbed 
lattices, and is one for MRJ. In Ref. Q we examined the possibility of using the surface 
term coefficient B as an order metric (increasing B indicated greater disorder). We did 
not anticipate the appearance of a further logarithmic term for the disordered packings. In 
this sense, the MRJ packings are markedly more disordered: they have macroscopic density 
fluctuations which are much larger than crystalline packings. Quantitative understanding 
of this aspect of disorder and its relation to density fluctuations remains a fascinating open 
problem. 
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Figure 2: (Color online) Structure factor for a 10^-particle packing {(f) = 0.642) and for a hard- 
sphere liquid near the freezing point {(p = 0.49), as obtained via two alternative numerical methods 
and also from the Percus-Yevick (PY) theory for the liquid DFT results are also shown 

over a larger range of K for a 10^-particle packing ((/> = 0.643). The left inset focuses on the 
range close to the origin, showing that while a parabola matches the liquid data reasonably well [ 
S{K) ~ 0.02 + 4 • \Q~^K^ according to PY theory, which is known to underestimate Sifi) ], it does 
not appear appropriate for the jammed packing for large-to-intermediate wavelengths [as obtained 
from Eq. ((21)]. The very linear behavior of the DFT data for the jammed packing in the range 
up to K/2ti < 0.4 is remarkable. The right inset shows c(r) convolved (smudged) with a narrow 
Gaussian [due to numerical truncation of >S'(A;)]. The peak at r = D is thus in fact (almost) a 
5-function. 
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Figure 3: (Color online) The variance cr^ as a function of the window radius for a 10^-particle 
packing. The uncertainty in the variance, as shown with error bars, is estimated to be of the order 
of cr^/\/M, where M = 10^ is the number of windows used for a given window. Also shown is 
the theoretically predicted dependence of the form AX^ + CX^ InX + BqX^, along with just the 
surface term BqX"^, which dominates the density fluctuations. 
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Figure 4: (Color online) The excess coordination for a 10^-particle packing, along with the best 
fit of the form ^ for the tail, and the estimated uncertainty. Statistical fluctuations overcome 
the actual correlations after x ~ 15. Averaging over nine samples only shrinks the magnitude of 
the fluctuations by three, without revealing additional information. The inset on the top uses a 
logarithmic scale, and the inset on the bottom shows the zeroth and first running moments along 
with their asymptotic values as estimated from the tail fit. Note that for the range of x shown 
explicit finite-size corrections are small (less than 5%). 
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